1/1 


AD-A18S  479 
UNCLASSIFIED 


SPECTRAL  MULTIGRID  METHODS  FOR  THE  SOLUTION  OF 
HOMOGENEOUS  TURBULENCE  PRO.  .  (U)  INSTITUTE  FOR  COMPUTER 
APPLICATIONS  IN  SCIENCE  AND  ENGINEERIN.  . 

6  ERLEBACHER  ET  AL.  JUL  87  I  CASE-87-45  F/Q  12/1 


NL 


AD-A189  475 


FILE.  COE) 

NASA  Contractor  Report  178341 
ICASE  REPORT  NO.  87-45 


ICASE 


SPECTRAL  MULTIGRID  METHODS  FOR  THE  SOLUTION 


OF  HOMOGENEOUS  TURBULENCE  PROBLEMS 


G.  Erlebacher 
T.  A.  Zang 
M.  Y.  Hussalni 


D71C 

^iiLECTEBl^ 
^  DEU  2  4  1987^^ 


Contract  No.  NAS1-I8107 
July  1987 


PisrArbUTioN' T,  ■ 

Approved  foi  p  i 

Distribution  d  •„ 


INSTITUTE  FOR  COMPUTER  APPLICATIONS  IN  SCIENCE  AND  ENGINEERING 
NASA  Langley  Research  Center,  Hampton,  Virginia  23665 

Operated  by  the  Universities  Space  Research  Association 


NASA 

National  Aeronautics  and 
Space  Administration 

Langley  Research  Cent* 

Hampton,  Virginia  23665 


87  12  16  2 : 


SPECTRAL  MULTIGRID  METHODS  FOR  THE 
SOLUTION  OF  HOMOGENEOUS  TURBULENCE 

PROBLEMS 


G.  Erlebacher  and  T.  A.  Zang 
NASA  Langley  Research  Center 

M.Y.  Hussaini* 

Institute  for  Computer  Applications  in  Science  and  Engineering 


New  three-dimensional  spectral  multigrid  algorithms  are  analyzed  and  implemented  to 
solve  the  variable  coefficient  Helmholtz  equation.  Periodicity  is  assumed  in  all  three  di¬ 
rections  which  leads  to  a  Fourier  collocation  representation.  Convergence  rates  are  the¬ 
oretically  predicted  and  confirmed  through  numerical  tests.  Residual  averaging  results 
in  a  spectral  radius  of  0.2  for  the  variable  coefficient  Poisson  equation.  In  general,  non¬ 
stationary  Richardson  must  be  used  for  the  Helmholtz  equation.  The  algorithms  developed 
are  applied  to  the  large-eddy  simulation  of  incompressible  isotropic  turbulence. 


*  Research  supported  by  the  National  Aeronautics  and  Space  Administration  under 
contract  No.  NASl-18107  while  resident  at  the  Institute  for  Computer  Applications  in 
Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23065. 


%,s 


iL'tLSLM  ut  «-*  t.1  k.t  u  u  I j  Li'L i'Li'U  t.i'Li't.i'L .u  Au  a 


I.  Introduction 

Spectral  multigrid  methods  [9,13,17,18]  combine  the  accuracy  of  spectral  discretizations 
with  the  efficiency  and  flexibility  of  multigrid  solution  techniques.  To  date  they  have 
been  implemented  in  an  exclusively  two-dimensional  setting,  with  applications  to  elliptic 
model  problems  [9,17,18]  and  to  compressible,  potential  flows  [13].  In  this  paper,  spec¬ 
tral  multigrid  methods  are  extended  to  three-dimensional  periodic  problems  and  applied 
to  the  large-eddy  simulation  of  turbulent  flow.  This  work  represents  one  realization  of 
the  prospective  large-scale  applications  of  spectral  multigrid  methods  that  were  discussed 
in  [19]. 

In  section  2,  the  concept  of  large-eddy  simulation  is  introduced  and  the  discretized 
Helmholtz  equations  are  formulated.  Section  3  discusses  several  multigrid  algorithms  suit¬ 
able  for  Poisson  and  Helmholtz  equations.  Numerical  results  on  the  model  problems  are 
discussed  in  section  4.  Finally,  the  multigrid  algorithms  developed  in  the  previous  sections 
are  incorporated  into  the  full  time-dependent  turbulence  simulation.  Numerical  results  are 
given  in  section  5. 


II.  Incompressible  homogeneous  turbulence 

Large-eddy  simulation  (LES)  models  the  small  spatial  scales  of  a  turbulent  flow  as  a  func¬ 
tion  of  the  large  scale  variables  [2,5,10,11,12,15].  A  spatial  filter  applied  to  the  velocities 
produces  the  large-scale  velocities  from  which  the  small  spatial  scales  have  been  removed. 
The  validity  of  LES  rests  on  the  assumption  that  the  small  scale  statistics  are  insensitive  to 
geometry  away  from  solid  boundaries.  Inasmuch  as  this  criterion  is  satisfied,  a  good  LES 
model  is  applicable  to  a  wide  variety  of  configurations.  Alternatively,  the  Reynolds  aver¬ 
aged  Navier-Stokes  equations  result  from  time  averaging  the  Navier-Stokes  equations  [14]. 
The  resulting  perturbation  velocities,  the  difference  between  the  true  velocities  and  the 
averaged  velocities,  become  the  velocity  fluctuations  in  time  and  contain  information  at 
all  spatial  scale  lengths.  Consequently,  Reynolds  averaged  turbulent  models  are  expected 
to  be  more  geometry  dependent  than  LES  models. 

In  non-dimensional  form,  the  conventional  Navier-Stokes  equations  are  given  by 

—  I  V  •  (t7v)  -  -  Vp  +  V  •  i/Vv  (1) 

VvO  (2) 

where  v  is  the  velocity  vector,  p  is  the  static  pressure,  and  i/,  the  kinematic  viscosity,  is 
assumed  to  be  constant. 

Any  flow  variable  7  can  be  spatially  filtered  in  the  following  manner: 

7(x)  J  1 7(x  z,  A)/(z)<i3e  (3) 
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where  G  is  a  filter  function,  A  is  the  computational  mesh  size,  and  D  is  the  domain  of 
the  fluid.  It  follows  that  Eq.  (3)  substantially  reduces  the  amplitude  of  the  high-frequency 
spatial  Fourier  components  of  any  flow  variable  7.  Consequently,  7  can  be  more  accurately 
termed  the  large-scale  part  of  7 . 

The  turbulent  fields  are  decomposed  into  their  large  and  small  scale  components  based 
on  the  prescription 

7  =  7+r  Hi 

where  7'  is  the  velocity  representative  of  the  small  spatial  scales.  The  direct  filtering  of 
the  momentum  equation  yields 

dv  _  — 

—  +  V-(t;t>)  =  -  Vp  +  V-i/V v  -I-  V-r  (5) 


Vu  =  0 


where 


Tu  =  -(fffcffj  -  -I-  V*uj  +  v[vk  +  (7) 

is  the  subgrid-scale  stress  tensor.  This  tensor  can  be  decomposed  into 

Lki  =  h)  (8) 

Cki  =  -(vfr  +  tfui)  (9) 

Rki  =  ~v'kvl  (10) 

which  are  respectively,  the  subgrid-scale  Leonard,  cross,  and  Reynolds  stresses  [6]. 

The  deviatoric,  i.e.  trace-free,  part  of  the  subgrid-scale  Reynolds  stress  tensor,  pR  ,  is 
approximated  by  the  Smagorinsky  model 

oRki  =  Ve  D$kl  (11) 

where  i/R  is  the  velocity  dependent  eddy  viscosity 

=2C*A2/4/j  (12) 

with 

—  1  dvk  dvt 

(13) 

th~  SmnSmn  (14) 

(i.e.,  S  is  the  Favre  filtered  rate  of  strain  tensor  while  I Ij  is  its  second  invariant)  and 
CR  is  the  Smagorinsky  constant  (the  Einstein  summation  convention  for  repeated  indices 
is  assumed.)  The  cross  and  Leonard  subgrid-scale  stresses  are  approximated  with  the 
Bardina  model  jl].  Together  with  the  subgrid-scale  Reynold  stresses,  one  obtains  the 
linear  combination  model  [l] 

Tkt  =  -~D{VkVl  -  VkVl)  +  I'E  pSkl)-  (15) 
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For  purposes  of  numerical  computation,  the  subgrid-scale  stresses  are  partitioned  into 
the  subgrid-scale  Reynolds  stress  and  the  remaining  terms.  The  latter  terms  contain  no 
derivatives  of  velocity,  and  are  therefore  treated  explicitly  along  with  the  advection  term. 

Substitution  of  the  subgrid-scale  stress  (15)  into  Eq.  (1)  transforms  the  momentum 
equation  into 

+  V-(«S)  =  -VP  +  V-(i/  +  uE)Vv  +  V-(L  +  C)  (16) 

where  the  isotropic  components  of  the  total  subgrid-scale  stress  have  been  lumped  together 
with  the  pressure  to  produce  a  new  pressure  variable,  P  ,  defined  by 

P  =  p  +  iLkk  +  iCkk  +  iRkk-  (17) 

The  subscript  I  indicates  that  only  the  trace  of  the  tensor  is  considered. 

For  the  isotropic  turbulence  problem,  equations  (16)  and  (6)  are  solved  in  a  cubic 
computational  domain,  periodic  in  all  three  spatial  directions.  Fourier  spectral  methods 
are  an  established  approach  to  this  problem  [10,12].  The  solution  is  obtained  in  two  steps. 
In  the  first  step,  the  convective  terms  and  the  Leonard  and  cross  subgrid-scale  stresses 
are  solved  explicitly  while  the  viscous  terms  are  treated  with  an  implicit  algorithm.  In 
the  second  step,  a  Poisson  equation  is  solved  for  the  pressure  to  insure  that  the  velocity 
field  remains  divergence-free.  For  convenience,  the  overbars  are  removed  hereafter  from 
the  primitive  variables,  and  it  is  understood  that  the  variables  refer  to  spatially  averaged 
quantities.  For  a  first-order  time  discretization,  the  first  step  thus  solves 

tT*  =  xT-  At[cD  x  t7+ V-(dL  +  DC)|"  +  A<V(i/ +  i/£)Vtf\  (18) 

Note  that  the  momentum  equation  is  used  in  rotation  form.  The  pressure  therefore  acquires 
the  additional  term  l/2|u|2.  As  a  result  of  the  first  step,  one  obtains  an  intermediate 
velocity  field  tT*  which  serves  as  initial  conditions  for  the  correction  stage 

0"+l  =  v"  —  At  VPM  1  (19) 

VtT+l  =  0  (20) 

For  spectral  collocation  algorithms,  a  direct  solution  of  Eq.  (18)  is  not  feasible  because 
the  matrices,  which  represent  the  diffusion  operators,  are  full.  The  alternative  which  is 
discussed  here  is  to  use  iterative  methods,  and  in  particular  spectral  multigrid  (SMG) 
methods,  for  the  solution.  The  second  step  in  the  splitting  algorithm  can  be  transformed 
into  a  Poisson  equation  with  constant  coefficients  which  is  solved  exactly  in  Fourier  spac'\ 
In  practice,  both  the  explicit  terms  are  solved  with  a  third  order  Runga-Kutta  algorithm, 
while  the  implicit  diffusion  terms  are  approximated  with  a  Crank-Nicholson  scheme.  The 
implicit  equations  are  solved  at  each  of  the  three  Runga-Kutta  stages. 
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III.  Spectral  Representation 

When  the  solution  to  a  numerical  problem  is  approximated  by  a  truncated  series  of  ap¬ 
propriate  global  basis  functions,  the  solution  is  said  to  have  a  spectral  representation. 
The  method  of  projecting  the  solution  onto  the  basis  function  space  determines  the  type 
of  spectral  approximation:  Galerkin,  tau  or  collocation.  Spectral  methods  are  explained 
thoroughly  in  [4,9]  and  a  summary  of  their  applications  to  fluid  dynamics  is  provided  in  [7], 
In  this  paper,  only  collocation  methods  are  considered,  since  they  are  better  suited  to  the 
solution  of  non-linear  and  variable  coefficient  problems.  Periodicity  further  restricts  us  to 
a  Fourier  representation  of  the  primitive  variables. 

Consider  the  three  dimensional  periodic  function  u(r)  on  the  domain  [0, 27rj3  and  its 
truncated  Fourier  representation 


,N„N„N.IA  = 


N,/2  N,/ 2  N./ 2 

E  E  E  (2i) 

k,  =  -N,/t+l  kt  =  -Nt/ 2+i  k.  =  -N./ 2+1 


The  superscripts  on  u,  henceforth  omitted,  refer  to  the  set  of  collocation  points 

=  (Zj,yk,zt)  =  (jhz,khv,lht),  0  <  j  <  Nz,0  <  k  <  Nv,0  <  l  <  Nz.  (22) 

where  h  =  (hz,hy,hz)  =  (2ir  /  Nz,2ir  /  Nv,2ir  /  Nz).  The  number  of  collocation  points  in  the 
x,y,z  directions  are  respectively  (Nz,  Ny,  Nz).  To  insure  spectral  accuracy,  u(r)  must  be  a 
C°°  function.  The  function  u,  evaluated  at  the  collocation  points  m,n,p,  and  the  Fourier 
coefficients  are  related  through  the  pair  of  discrete  Fourier  transforms 


N./2  N,/2 

E  E 

N./2 

(23) 

-N./2+1  kw  =  -N,/ 2+1  *. 

=  - N./2+1 

1  N’ 

N  N  N  ^ 
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n=  0 

P=  o 

(24) 

First  and  second  derivatives  of  u  are  simply  obtained  by  differentiating  Eq.  (21)  term 
by  term  and  evaluating  the  result  at  the  collocation  points.  For  example,  the  first  and 
second  derivatives  in  the  x-direction  are 

du  N./2  N,/2  N./2 

Q-  =  Y  Y  Y  ^k„kt,k.ei{klZm+k,Vn+k,tp)  (25) 

m,n,p  k,  =  -N,/ 2+1  k,  =  -Nt/ 2+1  k,  =  -N./2+ 1 


N./2  N,/2  N./2 

=  Y  Y  Y  (- kl)uk.,k„k,ei{k’:tm+k'Vn+k,tr) .  (26) 

m  n,p  k,  =  -N,/ 2+1  k,=  -Nt/ 2+1  *,  =  -/V,/2+l 
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Derivatives  in  the  y  and  z  directions  have  similar  expressions.  When  evaluating  first 
derivatives,  the  highest  mode  in  the  direction  the  derivative  (the  N/2  mode)  is  removed 
since  it  makes  a  purely  imaginary  contribution  to  the  first  derivative  at  the  collocation 
points. 

In  many  problems,  one  is  required  to  solve  large  systems  of  equations  on  fine  grids. 
However,  direct  methods  are  often  impractical  because  of  the  size  of  the  problem,  and 
standard  iterative  methods  have  very  slow  convergence  rates.  Typically,  the  high  frequency 
components  of  the  error  damp  out  quickly,  while  there  is  a  very  slow  decay  of  the  error 
on  the  larger  scale  lengths.  Such  relaxation  schemes  smooth  out  the  error  very  quickly. 
Multigrid  methods  accelerate  the  convergence  of  iterative  methods  by  recognizing  that  low 
frequency  errors  on  a  fine  grid  become  high  frequency  errors  on  a  coarse  grid.  Therefore, 
the  smoothed  residual  is  interpolated  onto  a  coarser  grid,  and  a  new  set  of  equations,  similar 
to  the  original  set,  is  solved.  The  coarsening  process  is  continued  until  a  sufficiently  coarse 
grid  is  reached  on  vhich  a  direct  solution  procedure  is  relatively  inexpensive.  From  the 
coarsest  grid  solution,  >  solution  on  the  next  finer  grid  is  obtained  by  prolongation  of  the 
coarse  grid  correction  onto  the  next  finer  grid,  optionally  followed  by  several  relaxation 
sweeps  to  eliminate  the  high  frequency  errors  introduced  by  the  interpolation  process.  In 
general,  therefore,  multigrid  algorithms  have  three  components:  a  restriction  operator  to 
transfer  residual  information  from  the  finer  to  coarser  grids,  a  prolongation  operator  to 
extend  a  coarse  grid  correction  to  the  next  finer  level,  and  a  smoothing  algorithm  whose 
objective  is  to  reduce  the  high  frequency  components  on  a  given  level.  There  exists  an 
extensive  literature  on  multigrid  algorithms.  Several  good  review  papers  appear  in  [8], 

Spectral  multigrid  distinguishes  itself  from  other  types  of  multigrid  approaches  in  the 
choice  of  the  interpolation  and  prolongation  operators.  In  the  problem  considered  here,  all 
functions  are  periodic.  This  leads  to  the  preferred  truncated  Fourier  series  representation. 
Following  [17],  interpolation  of  a  variable  from  a  fine  to  coarse  grid  consists  of  the  following 
steps.  Transform  the  variable  to  Fourier  space,  reject  the  highest  modes  not  resolvable  on 
the  coarse  grid,  and  transform  back  to  physical  space  on  the  coarse  grid.  Prolongation 
is  done  in  a  similarly  straightforward  manner.  After  transforming  the  variable  to  Fourier 
space,  additional  terms  are  added  to  the  Fourier  series  with  zero  coefficients.  The  newly 
defined  function  is  then  transformed  back  to  physical  space  on  the  fine  grid.  Contrary 
to  the  more  popular  interpolation  methods  used  in  the  finite-difference  context,  which 
always  introduce  high  frequency  components  into  the  solution,  the  spectral  interpolation 
just  described  is  exact  for  solutions  to  the  constant  coefficient  Helmholtz  equation. 

Fast  Fourier  transform  (FFT)  methods  permit  the  basic  interpolation  calculations  to 
be  performed  in  0(iV3log  N)  floating  point  operations,  when  the  number  of  nodes  in 
all  three  directions  is  equal  Nx  —  Nv  =  NM  =  N.  The  grid  transfer  operators  and  the 
residual  calculation  are  both  based  on  FFT’s,  the  former  to  interpolate  variables  between 
different  grids,  and  the  latter  to  perform  first  and  second  derivative  evaluations  at  the 
grid  points.  Therefore  the  overall  multigrid  scheme  has  an  operation  count  proportional 
to  0(N3  log  N). 


A.  Relaxation  Scheme 


The  use  of  FFT’s  to  calculate  the  residual  restricts  the  choice  of  relaxation  schemes  to  si¬ 
multaneous  relaxation  schemes  such  as  Jacobi  and  Richardson.  These  relaxation  schemes 
are  implemented  in  physical  space.  Consider  the  constant  coefficient  scalar  Poisson  equa¬ 
tion 

V2u  =  /(f).  (27) 

The  Richardson  scheme  is  one  of  the  simplest  smoothers.  Applied  to  Eq.  (27),  the  solution 
after  one  smoothing  step  becomes 

u  <—  u  —  wr,  (28) 

where  r  is  the  residual  /  — V2u  and  u>  is  the  relaxation  parameter.  Both  stationary  (fixed  w) 
and  non-stationary  (variables)  are  considered.  Equation  (28)  admits  a  Fourier  analysis.  If 
a  single  three-dimensional  Fourier  mode  ( jyk,l )  is  substituted  into  Eq.  (27),  the  smoothing 
rate,  jz,  becomes 

MW  =  |1  ~  "(hi  +  %  +  k])\  (29) 

where  w  is  the  relaxation  parameter.  In  the  context  of  multigrid  methods,  the  objective 
is  to  minimize  the  smoothing  rate  of  the  high  frequencies  seen  by  the  fine  grid,  and  not 
resolvable  on  the  coarser  ones.  Given  an  existing  grid,  the  next  level  of  coarsening  is 
obtained  by  defining  the  set  of  collocation  points  fljj,.  The  range  of  wavenumbers  over 
which  the  minimization  is  performed  is  the  difference  between  the  two  cubes  (in  wave 
number  space)  [0,  y]3  and  [0,  y]3.  Strictly  speaking,  because  both  the  mean  mode  and  the 
N/2  mode  have  been  filtered  out  of  the  right  hand  side,  /(f),  the  wave  numbers  considered 
for  the  minimization  should  actually  be  in  the  region  defined  by  the  difference  between  the 
cubes  [l,  y]3  and  [l,  y]3. 


A  straightforward  calculation  leads  to  an  optimal  smoothing  rate  of 


/z  = 


Nd-\ 

Nd+\ 


(30) 


for  the  Richardson  iteration  scheme,  where  Nd  is  the  number  of  spatial  dimensions.  When 
Nd  =  3,  /Z  =  11/13  «  .85.  It  is  obvious  from  Eq.  (30)  that  the  asymptotic  smoothing  rate 
increases  with  increasing  spatial  dimension.  For  example,  as  the  number  of  dimensions 
increases  from  1  to  3 ,  jz  increases  from  0.6  to  0.85.  For  3-D  problems,  the  optimal  relaxation 
parameter  is 


32 

137V2 


(31) 


Operators  that  are  spectrally  discretized  have  a  wider  spread  of  eigenvalues  than  their 
finite-difference  counterparts.  For  example,  the  optimal  Richardson  smoothing  rate  for  a 
second-order,  central  difference  discretization  of  the  constant  coefficient  Poisson  equation 
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In  contrast  to  the  spectral  smoothing  rates,  Pjrp  ranges  from  0.33  to  0.71  for  1,  2  and  3-D 
problems. 

Brandt,  Fulton  and  Taylor  [3j  applied  the  residual  averaging  technique  to  accelerate 
the  convergence  of  Richardson’s  smoothing  algorithm  for  two-dimensional  Fourier  repre¬ 
sentations.  The  extension  to  the  present  three-dimensional  problem  is  straightforward. 
The  smoothing  algorithm  now  satisfies 


umnp  —  umnp  jyj  ■^■rmnp 


where  (m,n,p)  is  the  grid  point  to  which  the  smoother  is  applied  and  A  is  the  averaging 
template.  A  is  the  three-dimensional  array 
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If  elements  of  A  are  denoted  by  Atjk,  the  three  arrays  in  the  above  expression  correspond 
(from  left  to  right)  to  k  =  1,2,3.  The  parameter  S  is  associated  with  the  8  corner  points 
of  the  cube  centered  at  (m,n,p)  about  which  the  averaging  is  being  performed.  In  con¬ 
ventional  notation,  the  averaging  template  A  applied  to  the  residual  at  (m,  n,p)  produces 
the  expression 

Arm  np  —  [a  +  0  S  J  rm+i,n+j,p+k-  (35) 

l<|+UI+l*l=l  |i|+l>|+|k|=2  l«l+|j|+|fc|=3 

Such  a  scheme  is  called  weighted  residual  averaging  (WRA).  A  Fourier  analysis  applied  to 
Eq.  (33)  yields 

fi(kt,kv,k,;a,0,i,6)  =  |1  -  ~{k\  +  k\  +  *’)[a  +  20{cosOx  +  cos 6V  +  cos 0t) 

+4q(cos  0V  cos  0M  +  cos  0X  cos  0X  +  cos  0Z  cos  0y) 
+8^(cos5Icos^1/cos^,)j|.  (36) 

where  ( 0X,0V,O ,)  is  a  shorthand  notation  for  %{kx,kv,kz).  The  solution  to  the  minimax 
problem 

p  =  min  max  n{0x,0v,Oz,a,0,i,6)  (37) 

yields  the  optimum  parameters  a,0,i,  and  6  as  well  as  the  smoothing  rate  p.  This  is  solved 
numerically.  The  angles  lie  in  the  region  formed  by  the  difference  between  the  two  cubes 
[0,7r/2]3  and  [0,  ?r/4]3.  To  demonstrate  the  importance  of  all  four  averaging  coefficients,  p 
was  calculated  by  successively  increasing  the  number  of  non-zero  coefficients.  The  results 
are  presented  in  table  1  where  a  comparison  is  made  with  the  2-D  results  of  Brandt,  Fulton 
and  Taylor  [3].  In  both  2-D  and  3-D,  p  decreases  substantially  with  the  help  of  residual 
averaging.  As  expected,  the  minimum  3-D  spectral  radius  is  higher  than  the  optimal  2-D 
value. 
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0.062  0.000  0.0 
0.101  0.145  0.0 

0.120  0.287  0.0083 
0.144  0.042  0.0185 


<5  I  3-D  n 


0.0  0.852 

0.0  0.608 

0.0  0.453 

0.0085  0.195 


Table  1:  Optimal  averaging  parameters  and  smoothing  rates  for  weighed  residual 
averaging  scheme. 


B.  Variable  Coefficient  Poisson  Equation 

The  analysis  of  the  previous  section  is  exact  for  the  constant  coefficient  Poisson  equation. 
More  generally,  one  wishes  to  solve 

V-a(r)  Vu  =  f{r)  (38) 

where  a(f)  is  a  strictly  positive  C°°  function.  Although  convergence  is  no  longer  theoret¬ 
ically  guaranteed,  good  results  are  obtained  if  the  residual  is  first  divided  by  a(r)  before 
averaging.  Alternatively,  one  can  use  Eq.  (33)  after  replacing  4n2/N 2  by  47r 2/a(r)N2. 
Both  approaches  yield  similar  results.  The  results  presented  herein  are  based  on  the  for¬ 
mer  method. 


C.  Helmholtz  Equation 

With  the  good  convergence  rates  achieved  for  the  Poisson  equation,  attention  is  now  fo¬ 
cussed  on  the  three-dimensional  Helmholtz  equation 

V-aVu  -  A u  =  f{r)  (39) 

where  a(r)  is  a  C°°  function  and  A  is  a  positive  constant.  The  iteration  scheme  and  its 
associated  convergence  rate  are  respectively 

u  <—  u  -  U)  (f(r)  -  V-aVu  +  Au)  (40) 

and 

M  =  |1  -  u(k2  +  A)|  (41) 

where,  as  previously,  the  number  of  grid  points  is  assumed  to  be  equal  in  all  directions 
( Nx  =  Nv  =  Nt  =  N).  The  optimum  smoothing  rate  is 

H/13  ,  . 

M  -  32  A 

+  13A*  a 


which  is  exact  for  constant  coefficient  a.  (It  is  assumed  that  the  kz  =  kv  =  kz  =  0  mode  is 
solved  for  exactly  on  the  coarsest  grid.)  Note  that  the  difference  in  smoothing  rates  of  the 
Poisson  and  Helmholtz  operators  decreases  with  increasing  grid  size. 


Experiments  were  also  performed  with  a  non-stationary  Richardson  smoothing  algo¬ 
rithm.  If  k  is  the  condition  number 

k  —  ■  (43) 


of  the  discrete  spectral  operator  V-aV  —  A,  the  jth  relaxation  parameter,  of  a  k- 
parameter  cycle  is 

2/ Amjn 


uk  = 

j 


(«  -  l)cos(-2j  ^7-—-)  +  {k+1) 


l,...,k 


(44) 


2k 


and  the  corresponding  smoothing  rate  is 


1  n 


(45) 


which  is  the  solution  to  a  standard  minimax  problem  [16].  The  range  of  frequencies  that 
are  preferentially  damped  are  in  the  domain  defined  by  the  difference  between  the  two 
cubes  [0,Amai]3  and  [0,  Am,„]3.  As  a  function  of  A  and  of  the  coefficient  a(r),  the  minimum 
and  maximum  eigenvalues  of  the  discrete  Helmholtz  operator  are 


A  mi  n  — 


aN 2 


+  A, 


.  _  3 aN2 

A  max  —  "t  A 


(46) 


16  4 

and  the  sequence  of  optimal  relaxation  parameters  in  the  non-stationary  Richardson  iter¬ 
ation  scheme  become 

32  1 


U’  aN2  +  16A 


(*  -  i)cos('"in.  +  (*  + 1) 


(47) 


2k 


The  number  of  terms  in  the  sequence  is  set  equal  to  the  number  of  smoothing  sweeps.  For 
stationary  Richardson,  j  =  1  and  Eqs.  (47)  reduces  to 

32 


u 


13  N2 


1  + 


32  A 


(48) 


13  N2  a 

while  Ji  is  given  by  Eq.  (42)  If  the  value  of  Amin  and  Amai  (Eq.  (46)  are  inserted  into  the 
condition  number  defined  by  Eq.  (43),  it  is  clear  that  the  convergence  rate  must  increase 
when  either  N ,  or  A  is  increased. 


Table  2  confirms  that  JI  decreases  with  increasing  A  and  j.  In  practice,  a  3-cycle  scheme 
is  sufficient  to  reduce  the  Li  norm  of  the  residual  by  a  factor  of  5.  An  alternate  formulation 
of  the  Richardson  method  is  obtained  when  the  Helmholtz  term  is  treated  implicitly,  i.c 

u-up(f(r)  -V-aVu) 


or  equivalently 


(50) 


u  <—  u  +  —  V-aVtt) 

where  uh  is  the  implicit  Helmholtz  relaxation  parameter 


uH  = 


Up 


1  +  A  up 

If  Up  is  the  optimum  relaxation  parameter  for  the  constant  coefficient  Poisson 
given  by 


32 


UJp  — 


1  O  \T2„  ’ 


(51) 

operator, 

f52l 


uh  is  identical  to  the  value  obtained  in  Eq.  (48).  Therefore  the  two  algorithms  are  identical 
for  all  a  >  0.  However  they  differ  from  one  another  for  non-stationary  Richardson.  Indeed, 
in  the  implicit  formulation,  one  chooses  the  uk  to  optimize  the  convergence  of  the  Poisson 
operator  which  leads  to  relaxation  parameters  independent  of  A.  The  A  dependence  is 
introduced  as  an  extra  positive  term  in  the  denominator  of  Eq.  (49).  This  is  in  contrast 
to  wk  given  by  Eq.  (47)  where  A  appears  in  the  coefficient  of  the  cosine  function.  Table  3 
illustrates  the  differences  between  the  two  approaches  when  non-stationary  Richardson  is 
used  with  cycles  of  varying  length.  The  two  methods  give  approximately  identical  smooth¬ 
ing  rates,  except  in  the  limit  of  large  A  where  the  implicit  method  gives  slightly  better 
performance.  From  a  practical  point  of  view,  it  is  cheaper  to  evaluate  the  acceleration 
parameters  for  the  implicit  scheme  because  a  factor  1  /a  can  be  factored  out  of  uk-  and 
combined  with  the  residual.  This  is  not  possible  for  the  explicit  formulas  which  depend  on 


r. 


k 

/i(A  =  0) 

/Z(A  =  10) 

/Z(A  =  50) 

1 

0.846 

0.826 

0.755 

2 

0.747 

0.720 

0.631 

3 

0.689 

0.661 

0.573 

4 

0.655 

0.628 

0.542 

5 

0.634 

0.607 

0.524 

6 

0.619 

0.593 

0.512 

Table  2:  Smoothing  rates  for  non-stationary  Richardson  iteration  applied  to  the 
Helmholtz  equation. 


D.  Implementation 


The  stationary  and  non-stationary  relaxation  schemes  were  implemented  in  a  simple  V- 
cycle  formulation  which  is  described  in  detail  in  [17,18]  for  the  two-dimensional  Poisson 


10 


323 

Grid  Size 

1  643  1 

1283 

0.652/0.651  0.654/0.654  0.655/0.655 

0.628/0.618  0.648/0.645  0.653/0.653 

0.542/0.511  0.621/0.610  0.646/0.643 


Table  3:  Comparison  of  convergence  rates  of  explicit  versus  implicit  non-stationary 
Richardson  iteration  algorithms. 


equation.  Results  are  based  on  a  fine  grid  of  32s  and  3  coarser  grids  of  163,  83  and  43.  A 
constant  number  of  relaxation  sweeps  on  each  grid  on  the  upwards  (7VU)  and  downward 
( Nd )  branches  of  the  V-cycle  was  found  to  provide  good  overall  smoothing  rates  for  all 
the  cases  that  were  considered.  An  effective  rule  of  thumb  is  to  decrease  the  residual  by 
an  order  of  magnitude  on  each  grid.  This  leads  to  Nd  =  2  for  WRA,  because  of  the  high 
smoothing  rates,  and  Nd  =  3  for  the  non-stationary  Richardson  (NSR)  schemes.  This 
corresponds  to  an  approximate  decrease  of  the  L2  norm  of  the  residual  on  the  order  of  0.2 
and  0.05  for  the  WRA  and  NSR  respectively  (per  Nd  fine  grid  relaxations).  In  all  tests, 
Nu  —  1.  This  is  because  the  variable  coefficient  a  introduces  high  frequencies  into  the 
residual  after  a  prolongation. 

All  the  computations  presented  were  performed  on  the  Cray  2.  Fast  Fourier  transforms 
are  coded  in  Fortran,  and  achieve  a  100  Mflop  rate  on  grids  on  643  and  above.  Timings 
may  fluctuate  by  10%  —  20%  for  identical  runs  due  to  system  load. 


E.  Results 


There  are  many  different  approaches  to  measuring  the  efficiency  of  a  multigrid  algorithm. 
Ultimately,  the  user  is  interested  in  the  total  CPU  clock  time  a  code  takes  to  complete 
execution.  However,  this  timing  is  strongly  computer  dependent,  and  on  a  given  system, 
the  programmers  skill  can  greatly  influence  the  results.  It  is  therefore  necessary  to  supple¬ 
ment  the  CPU  time  with  more  intrinsic  measures.  The  simplest  measure  is  the  smoothing 
rate  /Z  of  the  smoother  on  the  finest  grid.  Unfortunately,  Ji  does  not  take  into  account  the 
work  done  on  the  coarser  grids,  nor  does  it  account  for  the  time  spent  performing  grid 
transfers.  Alternate  criteria  are  required  to  take  this  work  into  account.  One  method  is 
to  calculate  the  ratio  of  L2  norms  of  the  residual  after  and  before  a  single  V-cycle,  and 
calculate  a  smoothing  rate  per  V-cycle  as 

-  fcliJ  (53) 


where  Nj  is  the  total  number  of  fine  grid  sweeps  during  one  V-cycle  (Nu  +  Nf)-  Although 
Jly  is  still  not  useful  as  an  accurate  measure  of  efficiency  since  it  does  not  take  residual 
transfers  into  account,  it  can  help  establish  whether  the  high  frequencies  seen  by  each  grid 
are  damped  at  the  same  rate.  If  p  and  pv  are  unequal  and  the  number  of  relaxations  is  the 
same  on  every  grid  (except  perhaps  the  coarsest),  the  frequency  content  of  the  error  vector 
is  unevenly  distributed,  and  the  smoothing  rates  will  be  different  on  each  grid.  Therefore, 
pv  is  useful  as  a  diagnostic  tool. 


A  better  measure  of  the  overall  algorithm  efficiency  is  obtained  from  the  number  of 
equivalent  fine  relaxation  sweeps  defined  as 


Neq  = 


CPU  time  per  V  —  cycle 
CPU  time  per  fine  grid  relaxation 


The  equivalent  convergence  rate 


Mr  — 


(54) 


(55) 


measures  the  decrease  in  the  residual  norm  per  fine  sweep  taking  the  total  multigrid  over¬ 
head  into  account.  Together  with  the  total  CPU  time  per  V-cycle,  the  performance  of  the 
algorithm  can  be  ascertained.  In  what  follows,  performance  is  measured  exclusively  by  p 
and  pT.  Processing  time  is  a  function  of  the  number  of  relaxation  sweeps  on  the  way  up 
and  down  the  V-cycle,  and  on  the  grid  size.  These  parameters  are  kept  constant  within  a 
given  table  except  in  table  6  where  the  effect  of  fine-grid  resolution  is  studied.  Therefore, 
the  decrease  in  the  residual  norm  after  a  fixed  number  of  multigrid  cycles  is  an  objective 
measure  of  the  algorithm’s  efficiency. 


As  explained  in  detail  in  [18],  the  N/2  Fourier  mode  must  be  filtered  out  of  the  residual 
every  time  it  is  computed.  In  one  dimension,  this  is  done  in  physical  space  by  projecting 
the  residual  function  onto  the  space  orthogonal  to  [  —  !)*,  j  =  1  ,NX.  In  higher  dimensions, 
a  sequence  of  1-D  filtering  operations  is  performed.  The  filtering  operation  consumes 
approximately  25%  of  the  residual  calculation  on  a  323  grid.  The  influence  of  vectorization 
is  clearly  seen  from  the  decrease  in  relative  time  spent  filtering  as  the  grid  size  increases. 
For  example,  as  N  increases  from  32  to  128,  the  percentage  of  time  spent  filtering,  measured 
with  respect  to  the  total  time  spent  calculating  the  residuals  (which  includes  the  filtering), 
decreases  from  25%  to  13%.  Equivalently,  the  percentage  of  time  spent  in  3-D  derivative 
evaluations  during  the  computation  of  the  residual  increases  from  67%  to  77%  as  the  grid 
size  increases  from  323  to  1283. 


When  solving  the  Poisson  equation,  the  mean  value  of  the  right-hand  side  of  the  equa¬ 
tion  must  be  filtered  out  to  insure  convergence  of  the  residual  towards  zero  [18].  This  is 
done  once  at  the  beginning  of  the  calculation. 


All  the  numerical  experiments  were  done  with  the  Helmholtz  equation  (with  A  set  to 
a  constant  value).  The  Poisson  equation  is  obtained  by  setting  A  to  zero.  The  coefficient 


In  all  cases  considered,  the  right  hand  side  of  the  Helmholtz  equation  is  calculated  to  insure 
that  the  exact  solution  is 

u,r(f)  =  sin(iVI7r  sin(i))  sin(iVy7r  sin(y))  sin(A^z7r  sin(a))  (57) 

The  factors  Nz,  Ny  and  Nt  are  included  to  insure  that  the  complete  spectrum  of  spatial 
scales  are  equally  represented  in  the  error  vector.  This  insures  that  /Z  =  JIV.  Such  a 
solution  however,  precludes  a  direct  comparison  of  the  computed  and  the  exact  solution 
because  uez  is  no  longer  well  represented  by  the  collection  of  Fourier  modes. 

Convergence  results  for  the  VVRA  are  presented  in  table  4.  The  smoothing  rates  ob¬ 
tained  for  the  constant  coefficient  Poisson  equation  are  lower  than  the  theoretical  predic¬ 
tions.  This  is  mainly  because  the  analytic  results  presented  in  Table  1  are  only  valid  in 
the  limit  N  — *  oo.  Although  a(r )  varies  by  a  factor  20  across  the  physical  domain  when 
c  =  1,  Ji  remains  close  to  the  optimal  value  of  0.2.  Experiments  have  shown  that  a  large 
degradation  of  JL  occurs  for  c  greater  than  1.5.  Seven  V-cycles  reduced  the  residual  10  to 
11  orders  of  magnitude.  Timings  indicate  that  the  number  of  equivalent  fine  relaxations  is 
3.45  on  a  323  grid.  This  is  larger  than  is  expected  from  the  sum  of  the  work  done  on  the 
sequence  of  grids  obtained  from  summing  the  geometric  series  1  +  (|)3  +  (§)3  +  (^j)3  -\ - 


Table  4:  Rates  of  convergence  for  Poisson  equation 


which  leads  to  1.13  equivalent  work  units.  The  discrepancy  between  the  theoretical  and 
numerical  results  are  due  to  the  time  spent  in  the  grid  transfers  which  are  not  included  in 
the  geometric  series,  and  the  inefficiency  of  the  Cray  2  program  on  the  coarser  grids. 

Large  eddy  simulations  of  turbulence  require  that  the  numerical  scheme  be  capable  of 
resolving  a  wide  range  of  spatial  scale  lengths.  Furthermore,  the  velocity  fields  have  a 
quasi-random  distribution  over  the  set  of  scale-lengths  that  survive  the  the  filtering.  In 
typical  large  eddy  simulations,  the  smallest  fluctuations  seen  by  the  LES  code  are  on  the 
order  of  4  fine  grid  cell  widths  (albeit  with  small  amplitudes).  It  is  therefore  important 
to  ascertain  the  influence  of  spatial  structure  in  a(r)  on  the  smoothing  rate  of  the  Poisson 
and  Helmholtz  operators.  To  this  effect,  the  definition  of  the  coefficient  c(r)  is  slightly 


modified  according  to 


a(r)  =  1  +  £e«*K*)+eo.(n.»)  +  eo.(n.«).  (58) 

With  e  =  0.5,  a  varies  from  1  to  about  10.  The  influence  of  ne  is  to  introduce  high  frequency 
content  into  the  coefficient  without  affecting  its  range.  In  other  words,  the  higher  nf,  the 
smaller  the  distance  over  which  a  assumes  its  maximum  variation.  Table  5  shows  that 
small  values  of  ne  do  not  adversely  affect  the  convergence  rate  of  the  iteration  scheme. 


n( 

) U 

/J-y 

1 

0.17 

0.17 

2 

0.18 

0.18 

4 

0.55 

0.55 

8 

0.58 

0.62 

Table  5:  Effect  of  high  frequency  content  of  a [r)  on  the  smoothing  rate 


However,  the  rate  of  convergence  rate  quickly  deteriorates  for  n£  >  4.  Note  that  when  the 
variation  of  a  becomes  too  rapid,  there  is  a  discrepancy  between  JI  and  Jiv.  This  probably 
indicates  that  there  is  a  frequency  imbalance  across  the  various  grid  levels,  due  to  the  grid 
transfer  operators  which  do  not  properly  interpolate  a(r)  onto  the  coarser  grids.  Similar 
experiments  performed  on  the  Helmholtz  equation  indicate  that  A  has  a  beneficial  effect 
on  the  smoothing  rate  as  ne  is  increased.  Of  course,  JI  is  higher  than  the  worst  Poisson 
result  since  residual  averaging  is  not  allowed. 

The  Helmholtz  equation  is  numerically  solved  for  several  values  of  e  and  A.  Non¬ 
stationary  Richardson  with  a  cycle  of  3  produces  the  results  displayed  in  table  6.  Pairs 
of  numbers  correspond  to  (e  =  0/e  =  0.5).  As  expected,  for  a  fixed  e,  JI  decreases  with 
increasing  A  due  to  the  reduction  of  the  condition  number 

^max 
^min 

— aN2  +  A 
16 _ 

-aN2  +  A  ’ 

4 

For  the  larger  values  of  A,  the  effect  of  non-constant  coefficient  a  is  more  severe.  Whereas 
JI  is  almost  unaffected  by  £  for  A  <  40,  the  differences  in  smoothing  rates  are  substantial 
for  A  >  100.  For  instance,  when  A  =  1000,  JI  is  approximately  3  times  larger  for  e  =  0.5 
than  for  e  =  0.0.  Similar  trends  occur  for  JiT.  However,  the  influence  of  £  on  Jlt  is  less 
dramatic  at  high  A. 


(60) 


(59) 


r10 

FH 


1 

10 

100 

1000 

0.68/0.68 

0.66/0.67 

0.49/0.62 

0.16/0.42 

0.87/0.88 

0.86/0.86 

0.77/0.83 

0.63/0.76 

3(~5)  / 1  (-4) 

l(-5)/2(-5) 

5(-9)/2(-6) 

l(-ll)/2(-9) 

Table  6:  Rates  of  convergence  for  Helmholtz  equation 


Optimal  multigrid  methods  produce  convergence  rates  independent  of  the  grid  size. 
This  is  tested  for  the  Helmholtz  equation  at  t  =  0.5  and  A  =  10.  Table  7  indicates  that 
both  Jt  and  pT  are  approximately  constant  over  fine  grid  sizes  that  range  between  32s  and 
128s.  In  all  cases,  the  coarsest  grid  level  is  43.  Computer  timings  for  the  calculations  are 
also  presented.  They  are  normalized  to  1  on  the  643  grid.  In  all  cases,  the  code  is  stopped 
after  a  fixed  number  of  V-cycles.  Increased  time  spent  on  the  1283  grid  relative  to  the  643 
agrees  quite  well  with  theoretical  predictions.  This  is  in  contrast  to  the  extra  70%  CPU 
time  spent  on  the  323  grid  than  allowed  for  by  the  0(N  log  N)  scaling.  Poor  vectorization 
on  this  grid  is  the  probable  cause  of  this  discrepancy. 


grid  size 

32s 

64s 

1283 

¥ 

0.67 

0.67 

0.68 

HT 

0.86 

0.85 

0.85 

IM| 

ITnT 

1.9(-5) 

1.5(-4) 

2.5(-5) 

CPU  time  (arbitrary  units) 

0.18 

1.0 

9.0 

0(N*  log  N)  (arbitrary  units) 

0.10 

1.0 

9.3 

Table  7:  Grid  independence  for  Helmholtz  equation,  e  =  0.5,  A  =  10 


F.  Large  Eddy  Simulation 


The  implicit  stage  of  the  LES  equations  requires  the  solution  to  the  set  of  three  scalar 
Helmholtz  equations  given  by  equation  (18).  Defining 

a  =  (6„ 

<  V  +  uE(v)  > 


W5 


m 

m  m  *•  m 


■Vs 


,V  V 


Sytf 

V 


/  V  v 


A  - » 

v/V 


m 


'V- 

>s'*W 


yy. 

,s*>: 


I 


aw 


(63) 


<  i/  +  uE(v)  >  At 

Equ.  (18)  reduces  to  the  three  scalar  Helmholtz  equations 

V  •  a(t/)Vi>(r)  -  Au  =  tf , 

where  xT  is  the  flow  velocity  after  the  explicit  step.  The  above  definitions  of  a(r)  and 
A  insure  that  <  a  >=  1,  which  allows  the  numerical  results  to  be  compared  against  the 
results  in  the  previous  sections. 


A  major  difficulty  expected  to  reduce  the  efficiency  of  the  multigrid  implementation, 
when  compared  with  the  theoretical  and  model  problem  results,  is  that  a(r)  is  now,  in 
effect,  a  random  function  of  the  spatial  coordinates.  Numerical  experiments  indicate  that 
the  multigrid  algorithm  fails  to  converge  for  A  below  a  certain  threshold.  In  the  current 
code,  this  threshold  is  A  ss  100.  Convergence  rates  and  timings  are  shown  in  table  8. 
Calculations  were  performed  on  a  323  grid. 


Although  equation  (63)  has  the  same  functional  form  as  the  model  equation,  the  effect 
of  a{f)  and  A  on  the  overall  multigrid  efficiency  relative  to  a  purely  explicit  scheme  is  not 
easy  to  determine.  One  reason  is  the  strong  dependence  of  these  parameters  on  the  filter 
width  A,  kinematic  viscosity  u,  and  Smagorinsky  constant  CE.  To  understand  better  the 
interelations  between  these  parameters  and  the  possible  gain  of  a  multigrid  strategy,  let 


n  __  Atdiff 
A  tadv 

where  A tadv  and  A tdlfj  are  respectively  the  maximum  time  steps  calculated  for  the  explicit 
advection  and  diffusive  terms.  The  accuracy  of  the  simulation  is  mostly  determined  by  the 
advection  terms.  Implicit  algorithms  are  consequently  most  favorable  when  the  diffusion 
time  step  is  much  smaller  than  Atadv  (i.e.  when  R  <<  1).  Stated  differently,  a  fully 
explicit  solver  is  cheaper  than  a  mixed  explicit/implicit  scheme  when  R  >>  1.  For  Fourier 
methods,  the  maximum  advection  Courant  number  of  the  third  order  Runga-Kutta  method 
is  0.6,  which  leads  to 

Ax 

A  tadv  =  0.6  —  .  (65) 

where  t;  is  a  representative  fluid  velocity.  (The  estimates  in  this  section  are  for  a  one¬ 
dimensional  problem.)  On  the  other  hand,  the  diffusion  time  limit  is 


Atdiff  =  0.25 


Ax2 


(66) 


v  +  uE 

In  this  discussion,  all  the  variables  are  assumed  to  be  constant.  Equations  (65)  and  (66) 
combine  into 

_  vAx  .  , 

R.  =  . - •  (67) 


f{v)CRn\  Ax1  4-  i/' 

Several  substitutions  have  been  made  to  arrive  at  Eq.  (67).  The  filter  width  A  has  been 
replaced  by  nA Ax  to  separate  the  computational  grid  size  from  the  actual  filter  width. 


Thus,  when  is  increased,  the  filtering  is  stronger,  and  more  high  frequencies  are  removed 
from  the  large-eddy  velocities.  The  velocity  dependence  of  the  Smagorinsky  model  is 
included  in  f(v)  whose  magnitude  is  a  slowly  decaying  function  of  nA. 


As  confirmed  in  table  8,  the  convergence  rate  improves  with  increasing  A.  However, 
according  to  Eq.  (62),  a  higher  A  is  the  result  of  a  decrease  in  Cj?  or  u.  Other  sources 
of  variation  are  not  considered  here.  Equation  (67)  therefore  clearly  indicates  that  a  higher 
A  reduces  the  gain  of  the  multigrid  scheme  over  the  purely  explicit  scheme.  This  is  borne 
out  by  table  8.  The  multigrid  code  performs  2  times  slower  than  the  explicit  scheme  when 
A  =  130  and  9  times  slower  when  A  =  1000  although  JI  has  dropped  from  0.58  to  0.12. 
Furthermore,  as  A  increases,  so  does  Z  which  explains  why  the  multigrid  code  performs  so 
poorly  (compared  to  the  explicit  scheme)  for  large  A.  The  variation  of  the  results  in  the 
table  8  is  not  uniformly  monotonic  as  a  function  of  A.  This  is  partially  because  timings 
are  a  function  of  the  load  on  the  Cray  2,  and  of  the  interaction  between  the  various 
independent  control  parameters.  For  example,  if  A  is  increased  through  a  smaller  value 
nA,  the  high  frequency  content  in  the  velocities  is  reduced  and  better  SMG  convergence 
rates  are  expected,  not  only  due  to  the  larger  Helmholtz  coefficient,  but  also  because  of 
the  smoother  velocity  fields.  On  the  other  hand,  only  the  former  cause  for  improvement 
is  remains  when  Cr  is  decreased. 

Finally,  a(r)  was  assumed  to  be  constant  in  the  above  analysis.  In  all  probability,  the 
oscillatory  nature  of  a(r)  will  also  strongly  influence  the  performance  of  the  SMG.  When 
A  drops  below  100,  the  implicit  scheme  fails  to  converge.  This  might  be  related  to  the 
highly  oscillatory  coefficients  which  appear  when  simulating  turbulent  flows. 


A 

impl.  code  vs.  expl.  code 
time/step  time/run 

<  100 

non-converged 

130 

0.58 

16 

2 

200 

0.75 

21 

2 

240 

0.26 

12 

2 

480 

0.21 

10 

3 

540 

0.20 

12 

4 

700 

0.20 

14 

7 

1000 

0.12 

12 

9 

Table  8:  Numerical  results  from  SMG  incorporated  into  the  multigrid  code.  Figures 
refer  to  5  complete  time  steps. 


More  efficient  Helmholtz  solvers  must  be  developed  before  spectral  multigrid  algorithms 


will  strongly  outperform  the  explicit  schemes  for  the  simulation  of  turbulent  flows.  The 
quasi-random  coefficients  in  the  LES  equations  also  call  for  new  spectral  interpolation 
procedures  to  improve  the  robustness  of  the  method. 

G.  Conclusion 

Three  dimensional  periodic  Poisson  and  Helmholtz  equations  have  been  solved  with  a 
3-D  spectral  multigrid  algorithm.  Convergence  rates  for  the  Poisson  problem  are  best 
when  weighted  residual  averaging  is  adopted.  The  spatially  dependent  coefficient  a(f)  was 
allowed  to  vary  by  more  than  an  order  of  magnitude  without  affecting  overall  convergence 
rates.  Although  weighted  residual  averaging  is  impractical  for  the  3-D  Helmholtz  equation, 
non-stationary  Richardson  is  a  viable  alternative  for  a  wide  range  of  a  and  A.  At  a  fixed 
amplitude  variation,  high  spatial  frequency  content  of  a  had  a  deleterious  effect  on  JI  for 
the  Poisson  equation.  This  is  unfortunate  since  for  turbulent  simulations  the  variables  are 
necessarily  osciwatory. 

The  algorithms  herein  were  successfully  incorporated  into  a  full  3-D,  non-stationary 
incompressible  LES  code.  It  was  found  that  in  the  range  of  parameters  examined,  the 
SMG  takes  at  least  twice  as  long  as  an  explicit  calculation.  This  is  part  due  to  the 
relatively  large  spread  of  eigenvalues  for  a  Fourier  collocation  algorithm.  Another  cause  is 
most  probably  related  to  the  strong  and  rapid  variations  of  a(r). 
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